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PERIODIC, SMALL-AMPLITUDE SOLUTIONS TO THE SPATIALLY 
UNIFORM PLASMA CONTINUITY EQUATIONS* 
by J. Reece Roth 
Lewis Research Center 

SUMMARY 

The coupled set of first-order nonlinear differential equations 
X = c 0 + c i X + c 2 y + c 3 X y + c 4 X 2 + c 5 y2 
y = Aq + A x X + A 2 y + Ag xy + A 4 x 2 + A 5 y 2 

are solved by an approximate method which gives y(t) in closed form for the particular 
case in which the variables x ant * y vary periodically in time, the coefficients CX and 
are real, and the peak-to-peak variation of x is of small amplitude relative to the mean value 
of x- The peak-to-peak amplitude of y, however, is not necessarily small compared with 
the mean value of y. These equations are a generalized form of Volterra's problem of two 
conflicting populations, in which x is the population density of the prey, and y is the popula- 
tion density of the predator, orvice versa. Similar equations may also be derived from the 
spatially uniform neutral and charged-particle continuity equations in the field of plasma 
physics. These equations then describe relatively low-frequency oscillations in the number 
density of electrons and neutral atoms in a partially ionized plasma. 

The functional dependence of y (t ) is given in closed form by Jacobian elliptic functions 
when the assumption of small amplitude is made (i. e. , Ax/xi« 1), and when the parameters, 
C^, A^, and the mean value of x, are such that periodic motion occurs. The amplitude, 
period, and waveform of the fluctuations in y are given in terms of the parameters, C^, A^, 
(which may be positive or negative), and x^- When periodic fluctuations are observed, it is 
possible to predict the amplitude, period, and waveform if C^, A.., andx^ are known (as is 
usually the case in plasma physics applications). It is also possible to infer relations among 
Cj and Aj if x^ and the period, amplitude, and waveform are known (as is usually the case 
in the predator -prey problem). When Cg = Ag = 0, theamplitude, period, and functional 
form of these Jacobian elliptic functions are derived for all cases in which periodic motion oc- 
curs. In the general case in which Cg * 0 and/or Ag # 0, it is shown that the closed -form 
solutions may also be expressed in terms of Jacobian elliptic functions, and the problem is left 
in a general form that may be solved for particular cases of interest. This procedure does not 
provide the waveform or amplitude of x(t), the quantity whose peak-to-peak amplitude is as- 
sumed to be small . 

*Some of the material contained herein has been submitted to the Journal of Mathematical 
Physics. 



INTRODUCTION 


The equations to be studied herein are 

x = c 0 + c ! X + C 2 y + c 3 xy + C 4 X 2 + C 5 y 2 (1) 

y = A 0 + A i x + a 2 y + a 3 X y + a 4 X 2 + a 5 y2 ( 2 ) 

(A.11 symbols are defined in appendix A . ) These equations are a generalized form of 
Volterra's classical problem of two conflicting populations (refs. 1 to 9). A less general 
form of these equations, comprising fewer terms on the right side, was originally 
studied by Lotka, in an attempt to formulate a mathematical theory of the behavior of two 
conflicting species in a state of nature (refs. 1 and 2), and of certain oscillatory chemical 
reactions (ref. 3). This work was somewhat later refined by Volterra (refs. 4 and 5), 

whose work is summarized by Davis (ref. 6). When a mathematical model of the 
predator -prey problem is formulated, and the terms of the differential equations are 
multiplied by the proper coefficients, the solutions of the equations are periodic (ref. 4), 
which is consistent with the observed periodic fluctuations in the population of predators 
and prey (ref. 5). 

The mathematical theory of coupled, first order nonlinear equations such as equa- 
tions (1) and (2) has been discussed extensively by Lotka, Volterra, Davis, and by 
M. Frommer (ref. 7). The nonlinear nature of these equations has precluded the obtain- 
ing of general solutions and has limited investigations of this problem to obtaining par- 
ticular solutions to simplified versions of equations (1) and (2), to studying the stability 
and existence of periodic solutions to the equations, or to obtaining a numerical solution 
to the equations for particular values of A ^ and C^. 

Most of the literature on this problem is concerned with determining the conditions 
under which the population of one or both species will remain bounded or under which 
periodic fluctuations of the populations will occur (refs. 6, 8, and 9). Little emphasis 
has been placed on determining the peak-to-peak amplitude, period, and waveform of the 
solutions to equations (1) and (2) because of the extreme difficulty of the exact, nonlinear 
problem. .Analytical solutions have been obtained for the linearized case in which the 
peak-to-peak amplitude of both x and y are small by comparison with their mean 
values (refs. 3 and 4). The waveforms are sinusoidal in this case. 

In the present analysis, analytical solutions to equations (1) and (2) were obtained 
under the assumption that the fluctuations of only one of the two variables are. small by 
comparison with its mean value. The class of solutions of interest in this analysis are 
those in which steady, periodic oscillation exists. Solutions which decay to a constant 
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value and those wherein either population increases without bound in violation of the 
small -amplitude approximation were not examined in any detail. 


APPLICATION OF EQUATIONS 

Physical Significance of Terms in Differential Equations for 
Classical Problem of Predator and Prey 

Consider the following pair of coupled, nonlinear first-order differential equations, 
in which x 1® the population density of the prey and y is the population density of pred- 
ators: 


i = C 0 + C lX + C 2 y + C 3 xy + C 4 x 2 + C 5 y 2 (1) 

y = A o + A 1 x + a 2 y + A 3 x y + A 4 X 2 + A 5 y2 (2) 

Equations (1) and (2) are a generalized form of Vol terra’s problem, in which the 
terms on the right side of each equation have the following physical significance: The 
terms Aq and Cq represent, respectively, the rate at which the predator or prey 
migrates across the boundary of the region of interest, when this rate is constant. The 
terms containing Ag and C ^ represent the e-folding growth rate at which the predator 
or prey would increase in the absence of the other species, in the absence of a constant 
rate of migration, in the absence of competition among the members of each species, but 
in the presence of an unlimited food supply for both species. The terms containing A^ 
and Cg represent the rate at which the predator and prey interact, and this rate mea- 
sures the rate at which the population of prey is depleted by the total number of predators. 
It is further assumed for mathematical convenience that the population of predators is 
augmented by the amount A^xY immediately upon the act of predation, which diminished 
the population of prey by an amount equal to C^xY- 

The terms containing C 4 and Ag represent the effects of interaction among the 
members of each species. The terms containing C 2 and A^ represent the beneficial 
or harmful effects on one species of the mere presence of the other, effects that are dif- 
ferent from those of predation and might, for example, consist of a form of symbiosis. 
Finally, the terms containing A^ and Cg represent the effects on the population of 
one species of competition among the members of the other species. This might come 
about, for example, if the prey fed in part on the carcasses of the predator. 
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One of the major problems in setting up a mathematical model of the predator-prey 
interaction is that there is often reason to suppose that the coefficients Cj and A i are 
themselves time dependent, varying with a seasonal periodicity (refs. 2 and 4). This 
variability greatly complicates the problem and makes the obtaining of closed-form solu- 
tions much more difficult. It is assumed in this analysis that the coefficients and 
are constant in time. It is also possible that more complex interactions between and 
among the two populations may result in a less simple interpretation of the individual 
coefficients, even when the total process is accurately represented by equations of the 
given form. 


Physical Significance of Terms in Differential Equations 
in Plasma Physics Applications 

Restricted forms of equations (1) and (2) have appeared in the literature of plasma 
physics (ref. 10), in which one equation was the energy equation and the other a contin- 
uity equation. Reference 10, however, contains no closed-form solution to these equa- 
tions. Apparently it has not been realized that if x is identified with the neutral density 
and y with the ion density, the plasma continuity equations themselves can be written in 
the form of equations (1) and (2), and the periodic solutions identified with a previously 
unrecognized plasma instability (ref. 11). 

Equations (1) and (2) can be derived in a plasma physics context by writing the con- 
tinuity equations for each of the three components (ions, electrons, and neutrals) of a 
partially ionized gas and assuming that the spatial variations are small over the region of 
interest. One of the three equations may be eliminated by assuming that the Debye dis- 
tance is small compared with the apparatus dimensions, so that the ion and electron 
densities are equal at each point. There will remain (in the absence of three -body pro- 
cesses) two equations similar in form to equations (1) and (2). The variable x might 
represent the number density of neutral particles and y the number density of electrons, 
or vice-versa. 

If the former is true, the terms in equations (1) and (2) containing Cq and Aq 
represent the rate at which neutrals and charged particles are injected or extracted from 
the volume of interest. The parameters C ^ and A^ represent the e-folding rate at 
which the density of neutrals and ions would change in the absence of other processes as 
a result, for example, of effusion across the boundary of the containment volume. The 
term containing Cg accounts for the decrease of neutral density due to ion-neutral or 
electron-neutral ionizing collisions. The term containing Ag gives the rate of increase 
of the ion density due to ion-neutral or electron-neutral ionization processes; this term 
also accounts for the scattering by neutral atoms or molecules of ions and/or electrons 
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into the escape cone of magnetic mirror -type confinement geometries. The terms con- 
taining and A ^ account for the change in ion or neutral number density due to 
neutral -neutral ionization processes, one neutral of which may be in an excited state. 

The term containing Cg gives the rate of increase of the neutral density due to ion- 
electron recombination processes. The term containing the parameter Ag accounts for 
the decrease in ion density due to ion-electron recombination processes, and also for the 
scattering into the loss cone by electron-electron, electron-ion, or ion-ion Coulomb inter- 
actions. The terms containing Cg and A^ represent the rate at which the number den- 
sity of one species changes because of the mere presence of the other. These latter two 
terms therefore do not appear to have physical significance in a plasma physics context, 
but they are included in the mathematical analysis for the sake of generality. The coef- 
ficients Cj and A- will usually be constant in plasma physics applications if the mean 
particle energy is independent of time. 

The mathematical analysis developed herein is sufficiently general so that the peak- 
to-peak amplitudes of x(t) and y(t) may or may not have the same magnitude. However, 
if x(t) refers to the density of neutral atoms or molecules in a partially ionized plasma 
and y(t) to the ion density, it is plausible that in many plasma physics applications the 
peak-to-peak amplitude of their fluctuations are about equal. This equality arises from 
the fact that ionizing collisions result in a one-to-one correspondence between neutrals 
lost and ions gained. The assumption that Ax/x^ « 1, therefore, also implies that 
Ay « Xj • The peak-to-peak fluctuation in ion density must then be small by comparison 
with the neutral density. This state of affairs is assured if the percentage of ionization 
is low at all times. 

The interpretation of x as die ion density and y as the neutral density is also pos- 
sible. In this case it would still be true that Ax « Ay, and the assumption that 
Ax/Xi « 1 would imply that Ay « x^ or that the change in neutral density is small by 
comparison with the total charged-particle density. The mathematical theory developed 
herein can therefore be applied to the case of a spatially uniform plasma, the percentage 
of ionization of which remains either low or high at all times. 

A detailed and extensive investigation of particular plasma physics applications of 
equations (1) and (2) is beyond the scope of this analysis. It is intended herein only to 
lay the mathematical groundwork required to apply equations (1) and (2) to plasma physics. 


ANALYSIS 

Objectives and Assumptions 


For an oscillatory system, a closed-form solution can show the relative importance 
of the various terms, as well as predict the effect of changes in the system on the char- 



acter of the oscillation. It is an easy matter to put equations (1) and (2) on a computer 
and obtain numerical solutions for a particular set of and A^. However, it is most 
desirable to have a closed-form solution for x(t) and y(t) in terms of the coefficients 
C i and A. of the differential equations, in order to compare experiment with the pre- 
dictions of the mathematical model. With such a solution the trend of the period and 
amplitude of the oscillation can be predicted as a function of experimentally determined 
quantities, which are contained in the constant coefficients and A^. 

The nonlinear nature of equations (1) and (2) defeated an attempt to find an unrestrict- 
ed closed-form solution to this set of equations. The most obvious approach to obtaining 
a closed-form solution is to linearize these equations through a ’’double-perturbation” 
analysis, in which both x and y are assumed to have constant mean values, Xj and 
yQQ. The time-varying portions of x and y ^(t) and are assumed to be much 

smaller than their respective mean values, so that x^ » and Yqo » y oi(^)- This 

method was used by Lotka (ref. 3) and by Volterra (ref. 5) in their analyses of a less gen- 
eral form of equations (1) and (2), and yields exponential and sinusoidal solutions. 

This double -perturbation approach, which is illustrated in appendix C, has at least 
two disadvantages. The first is that, in many interesting physical situations, the peak- 
to-peak amplitude of either x or y may not be small compared with the respective 
mean value. The second disadvantage of the double-perturbation approach is that the 
amplitudes X 2 ^) and yQ^(t) must be regarded as given small quantities, which cannot 
be obtained as functions of the coefficients Cj and A^. While the double -perturbation 
approach is perhaps more familiar and mathematically simple, the two disadvantages 
cited are deemed sufficient reason to study the more general problem discussed in the 
next paragraph. 

The small -amplitude approach used in this analysis starts with the assumption that 
the peak -to -peak amplitude of only one of the two variables is small by comparison with 
its mean value. Removal of this constraint from the second variable is a desirable gen- 
eralization; it permits the amplitude of the second variable to be calculated as a function 
of the coefficients and A^. This greater generality is achieved at some cost in 
mathematical complexity, since the resulting perturbation equations must be solved in 
terms of Jacobian elliptic functions rather than the real or imaginary exponential functions 
that suffice for the double -perturbation approach. In addition, the amplitude and wave- 
form of the first variable (whose peak -to -peak amplitude is assumed to be small) cannot, 
in general, be obtained in closed form. The period of both variables is, of course, the 
same. To recapitulate, the double -perturbation approach provides the waveform of both 
variables but the amplitude of neither; while the small -amplitude approach considered 
herein provides both the amplitude and the waveform of the variable with arbitrary peak- 
to-peak amplitude, but neither the amplitude nor the waveform of the variable whose 
peak -to -peak amplitude is assumed to be small. 
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An approximate closed-form solution is obtained herein by assuming that the peak- 
to-peak amplitude of the fluctuations in x is small compared with the mean value of x 


6 S M«1 

*1 


( 3 ) 


as is illustrated in figure 1. 




(b) y Q minimum. 

Figure 1. - Schematic illustrations of periodic solutions. 


The period of the waveform is defined as the time elapsed between two successive 
maxima of the functions x or Y- In the following analysis, it is assumed in all cases 
that the function y is periodic in time and that the boundary conditions at time t = 0 are 
given by equation (4). 


y(t = 0) = y 0 y(t = 0) = 0 


(4) 
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It is further assumed that Yq is positive definite and may be either the maximum value 
(fig. 1(a)) or the minimum value (fig. 1(b)) of y (t). 

No assumption was made about the relative magnitude of Ay and y, so that the as- 
sumption Ay/xi « 1 could also cover the case in which the relative amplitude Ay is 
either comparable to or small by comparison with the maximum value of y. In addition, 
the mathematical development is sufficiently general so that x^ may be either larger or 
smaller than the mean value of y. 

It appears that these assumptions, while restricting the range of application of the 
solutions, still cover an interesting and significant range of phenomena and should per- 
mit a contact to be made between theory and experiment in physical situations in which 
equations (1) and (2) apply. There are, however, certain limitations in the applicability 
of equations (1) and (2). One important limitation lies in the assumption that the coeffi- 
cients C i and Aj are not functions of time. In the predator-prey problem, not only do 
many processes have a seasonal dependence, but there is also a time lag between the 
appearance of a new member of the population and its participation in reproduction. 
Another limitation lies in the assumption that all the relevant processes can be included 
in a limited number of independent terms of no higher than quadratic degree. Because of 
this limitation, three-body recombination processes must be excluded from consideration 
in plasma physics applications. 


Development of Small-Amplitude Approximation 


The development was kept as general as possible, and the small -amplitude approxi- 
mation introduced at the last possible stage in the analysis. This procedure is both 
necessary and desirable, since information regarding the behavior of the variables is 
lost by introducing an approximation and this information does not appear in subsequent 
steps of the analysis. Premature introduction of the small -amplitude approximation 
therefore may result in obtaining a trivial or specialized solution. 

The basic equations of concern herein have already been given. 


x = C 0 + C 1 X + C 2 y + C 3 xy + C 4 x 2 
y = A Q + Aj x + A 2 y + A 3 xy + A 4 x 2 


+ c 5 y 2 

(i) 

+ a 5 y2 

(2) 


Immediate introduction of the small -amplitude approximation, by setting x(t) = x^ 

+ Ax(t) and then obtaining a solution for y(t) by ignoring terms containing Ay(t), would 
be premature since it would merely decouple equation (2) from the variation of x(t) and 
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lead to a family of trivial solutions for y(t). Of interest is the more general case, in 
which neither x nor y is zero, in which the waveform of y(t) depends on x(t) and on 
the coefficients (which in part determine the magnitude and time behavior of Ax(t)), 
as well as on the coefficients in equation (2). If equation (2) is differentiated with 
respect to time, 

V = A x x + A 2 y + A 3 xy + A 3 xy + 2A 4 xx + 2A & yy (5) 

The first-order time derivatives in equation (5) may be replaced by substituting 
equations (1) and (2) for them, to obtain 

y = [fiC 0 + A 2^o) + ^1 C 1 + A 1 A 2 + A 0 A 3 + 2A 4 C o) x 

+ (^1 C 4 + A 2 A 4 + A 1 A 3 + 2A 4 C l) x2 + ^3 A 4 + 2A 4 C 4) x3 ] 

+ [^1 C 2 + A 2 + A 3 C 0 + 2A 5 A o) + ( A 1 C 3 + 2A 2 A 3 + A 3 C 1 

+ 2A 4 C 2 + 2A 5 A l) x + ( A 3 C 4 + A 3 + 2A 4 C 3 + 2A 4 A 5) x2 ] y 

+ [( A 1 C 5 + 3A 2 A 5 + A 3 C 2) + ( A 3 C 3 + 3A 3 A 5 + 2A 4 C 5) x ] y2 

+ [ A s c s + ^fjy 3 (6) 

in which the coefficients of the powers of y(t) are grouped in brackets and the coefficients 
of various powers of x(t) in parentheses. 

Equation (6) is still completely general, but it is even more difficult of exact solution 
than the original equations (1) and (2). The only two time-varying quantities appearing in 
equation (6) are y(t) and x(t)* When the variable x(t) is written as x(t) = X^ + Ax(t), 
the terms within each bracket can be separated into two parts, one containing constants 
and the other depending on Ax(t). When the sum of the constant terms in each bracket 
is much larger than the sum of the time-varying terms containing Ax(t) (bearing in 
mind that Ax(t) « xp, these small terms may be considered negligible and omitted. 

This requirement is not automatically satisfied by the small -amplitude assumption for 
all possible cases but usually follows for particular numerical values of Aj, Cj, and x^- 
It is also possible that some of the expressions in parentheses or even those in brackets 
are negligibly small by comparison with other expressions of the same type. If so, their 
inclusion or exclusion is immaterial. 
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The variation of y(t) is expected to be dominated by the large, constant terms in 
equation (6) and only slightly perturbed by the small time-varying part containing Ax(t). 

It should therefore be possible to obtain an approximation to y(t) by neglecting the small 
time-varying terms containing Ax(t). Equation (6) thus may be approximated by an equa- 
tion in the single variable y(t) by introducing in it the small -amplitude approximation 
given by equation (3). 

If it is assumed, therefore, that the peak-to-peak variation of x(t) is sufficiently 
small that x(t) i n equation (6) may be approximated by its mean value x 

X = Xi + Ax « Xj (?) 

Equation (6) may then be written in terms of the constants A i? Cj, and x^ as 

y » [(AjC 0 + A 2 A 0 ) + (A 1 C 1 + A x A 2 + ApAg + 2A 4 C Q )x 1 

+ (a x C 4 + A 2 A 4 + A : A 3 + 2A 4 Cj)x5 + (a 3 A 4 + 2A 4 C 4 )xiJ 

+ [( A 1 C 2 + A 2 + A 3 C 0 + 2A 5 A o) + ( A 1 C 3 + 2A 2 A 3 + A 3 C 1 
+ 2A 4 C 2 + 2AgAj) Xl + (A 3 C 4 + A 2 + 2A 4 C 3 + 2A 4 A 5 ) X 3 ]y 

+ [( A 1 C 5 + 3A 2 A 5 + A 3 C 2,) + ( A 3 C 3 + 3A 3 A 5 + 2A 4 C 5) x l] y2 

+ [ A 3 c 5 + 2 4]y 3 ( 8 ) 

Since the terms in brackets are constant for a given problem, equation (8) can be written 
in the simplified form 

y = * 0 + l i y + h y2 + h y3 ( Q ) 

where the constants l . are defined by 

1 0 = [( A 1 C 0 + A 2 A o) + X 1^1 C 1 + A 1 A 2 + A 0 A 3 + 2A 4 C o) 

+ X i i C 4 + A 2 A 4 + AjAg + 2A^C^j + + 2A 4 C 4 ^J (10) 
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1 1 s [( A 1 C 2 + A 2 + A 3 C 0 + 2A 5 A o) + x l( A l C 3 + 2A 2 A 3 

+ AgCj^ + 2A 4 C 2 + 2AgA^ + x 2 1 (A 3 C 4 + A 2 + 2A 4 Cg + 2A 4 A^j (11) 

l 2~ [( A 1 C 5 + 3A 2 A 5 + A 3 C 2) + X 1^3 C 3 + 3A 3 A 5 + 2A 4 C s)] 

Z 3 = [( A 3 C 5 + 2A 5)] < 13 > 

The small -amplitude approximation has therefore reduced the general problem of finding 
the closed-form solutions of xOO and y(t) to the assumption that x(t) is approximately 
equal to xp y(t) is given by the solution to equation (9). The merit of the procedure 
used in the preceding calculations is that equation (9) happens to be one of the few non- 
linear differential equations whose solutions are available in closed form. Equation (9) 
is recognizable as the differential equation whose periodic solutions are given in terms of 
Jacobian elliptic functions (ref. 6, p. 129). 

The analytical procedure following equation (5) may be regarded in an equivalent but 
conceptually different light. In equation (5) xOO can be approximated by its mean value 
X^. Equations (1) and (2), with x(t) approximated by can then be regarded as 
approximations to x and y, which may also be substituted in equation (5). These 
approximations to x(t) (x and y in eq. (5)) result in an approximate equation for y 
that is algebraically identical to equation (9). 

Equation (9) may be integrated by multiplying both sides by y(t) and eliminating dt 
to obtain (ref. 12) 


d(y 2 ) = 2^ o + 1 1 y + 1 2 y2 + 1 3 y 3 ) d y ( 14 ) 

When this expression is integrated and the assumed boundary conditions given inequation (4) 
(y = yQ and y = 0 at t = 0) are applied, equation (14) becomes 

y 2 - 21 0 (y - *o) + 1 i(y 2 - 4) + f ■ 1 2 (y 3 - *o ) + ~ 1 - 4) < 15 > 

This equation is of a relatively well-known type (ref. 6, p. 129), whose solutions are either 
exponential or sinusoidal when 1 2 = 1 3 = 0. The solutions are given by Jacobian elliptic 
functions when 1 2 * 0 and/or 1 3 * 0, provided that the constant coefficients l . are real 
and are nondegenerate such that periodic motion occurs. Even the less interesting non- 
periodic solutions have implications for the long-term behavior of the system when the 
f. are such that it is not oscillating periodically. Therefore, all solutions to 
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equation (15) for real l j were investigated. 

In the following analysis, it is convenient to shift the origin of the variable y so that 


z = y -y 0 
z = y 


and 


Z(t = 0) = 0 
Z(t = 0) = 0 


(16) 

(17) 


(18) 


It should be noted that Z can be positive or negative, depending on whether y Q is the 
minimum or maximum value of y. The analysis is concerned with finding the total peak-to- 
peak amplitude of the fluctuation of y, 


z = v - y • 
max J max J min 

the functional form of Z(t), and the period T of the oscillation, where 


(19) 


Z (t = 0) = Z (t = T) (20) 

It is assumed that y(t) is real and positive but may be zero at one point on its waveform. 


Solution of the Case 1 3 f 0 

The most general solution to equation (15) occurs when 0. Referring to equa- 

tions (13), (1), and (2) reveals that this term is nonzero only if the y 2 term is present 
in either or both equations (1) and (2). 

If x(t) is regarded as the population density of the prey and y(t) as the population 
density of the predator, then l g is nonzero only when competition among the predators 
significantly affects their own population or directly affects that of the prey. In plasma 
physics applications, if x(t) is the density of neutral particles and y(t) is the charge 
density, is nonzero only if direct ion-electron recombination is a significant process 
or if significant numbers of collisions of two charged particles result in the loss of one 
of them, as by knocking one into an escape cone in a magnetic mirror confinement 
geometry. 
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Removing a factor of y - y Q from equation (15) leaves 

y 2 = (y - y 0 )[ 2 * o + 1 i( y + yo) + \ l 2 ( y2 + yy o + y o) + \ l s ( y3 + y2y o + yy o + y o)] < 21) 

Using equations (16) and (17) to convert equation (21) into an equation in the variable Z 
gives 


Z 2 = Z 


^(^0 + ^l y O + ^2 y O + ^3 y O + ^3 y o) + (^l + 2 ^ 9yn + 2 ^° y ”^ 


2 3 0 T 3 y 0> 

y o f + \ 


+ l -'2 + 2 ! 3 y 0 ) z 2 + ^ 3 z3 


Defining the auxiliary constants 


4 (^ Z 2 + Z 3 y o ) 


( 22 ) 


(23) 


. 2 ( t 1 + 2g 2 y Q +3g 3 y o) 


(24) 


and 


h = 4 (* 0 + 1 1 y 0 + ^2 y Q + ^3 y o) 


(25) 


makes it possible to rewrite equation (22) in the much simpler form 


Z 2 = ii 3 z(z 3 + jZ 2 + fZ + h) 
2 


(26) 


When the boundary conditions of equation (18) are used, this equation may be rewritten 
in integral form as 




1/2 


/ 

Jo 


dZ 


z ( : 


Z u + jZ + fZ + h 


■)] 


1/2 


(27) 
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This integral equation may be written in terms of the four roots of the denominator as 



where = 0 in all cases, and the roots, Z 2 , Zg, and Z 4 , are functions of the constant 
parameters, j, f, and h, given by the three roots of the cubic equation 

Z 3 + jZ 2 +fZ + h = 0 (29) 

The integral of equation (28) is immediately recognizable as being of a general type 
that gives rise to Jacobian elliptic functions (see, e.g., section 260 of ref. 13). 

The particular elliptic function involved and the expressions for the maximum amplitude 
of the fluctuations in y, the period of oscillation, and the elliptic modulus k all depend 
on the sign of l g and on the relation of the range of integration to the values of the four 
roots, Z^, Z 2 , Zg, and Z^. An exhibition of all possible solutions to equation (28) 
would be extremely tedious and is not attempted. In a specific instance in which l g # 0, 
however, it should not be difficult to determine the roots of equation (29) once the sign 
and magnitude of j, f, and h are known. It is then a relatively easy matter to integrate 
equation (28) (for which ref. 13 is most helpful) and to determine whether or not finite, 
periodic solutions exist and their specific form. 


Solutions of the Case ^ = 0 

The next most complex case is that in which 1 2 # 0 but l g = 0. The values of l q 
and l j are initially assumed to be arbitrary. In this case, equation (15) may be written 

y 2 - 21 „(y - y 0 ) ♦ i ,(y 2 - yj>) + 1 ‘ 2 (y 3 - yj) ( 30 > 

Factoring out y - Yq yields 

y 2 = (y - yq)[ 2z o + 1 i( y + yo) + \ l 2k 2 + yy o + y o)] ( 31 > 

Substituting equations (16) and (17) into equation (31) gives 

- z [f ' 2 z2 + (' i + 2> + 2 k o + yrf i + 2)] < 32 > 
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Defining two additional constants 


b s 


_ 3 ( z i + 2 y o Z 2) 


21 , 


c = 


_ 3 ( z o + y o z 1 + y o z 2) 


enables equation (32) to be written in the simplified form 


Z 2 =^ 2 z(z 2 +bZ + c) 


Equation (35) may also be written in terms of the roots of the right side as 


Z 2 =|! 2 (z-Zi)(z- z 2 )(z-Z 3 ) 


where Zg > Zg and 


and 



(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 


In order that the amplitude of Z have real values, the signs of l g, Z, and the roots in 
equation (36) must always combine to give a positive definite value to the right side of 
equation (35). The solution to the differential equation of equation (36) under the boundary 
conditions of equation (18) is given by the integral 
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The integral, or the right side of equation (40), is immediately recognizable as an 
elliptic integral of a type tabulated by Byrd and Friedman (ref. 13). Unfortunately, the 
functional form of the solution depends critically on the relation of the range of integra- 
tion to the roots, Zp Zg, and Zy A variety of solutions to equation (40) are presented 
in appendix B, and their ranges of validity and the nature of their solutions are indicated 
on table I. For the sake of concreteness, a bounded, periodic solution to equation (40) 
(case N of appendix B) is presented in the following paragraph. 

When both of the roots given by equations (38) and (39) are positive and the range of 
Z is confined to Zg ^ Z < 0, the solution to the basic integral of equation (40) is given 
by formula 233. 00 of reference 13. The argument of the elliptic function is given by 

t / 2 i ,\ 1/2 

U '*’ t)m ivn <41) 


where 


g 



The amplitude of Z has the functional dependence 


(42) 


Z = Z 2 SN 2 Ul (k) = - - 
2 


b + (b 2 - 4c) 


1/*1 


SN' 


u t (k) = y - y 0 


(43) 


The total peak-to-peak amplitude of Z is then 


Z = Z- = - ± 
max 2 g 


b + (b 2 - 4c) 


\l/2“ 


(44) 


and the modulus of the elliptic function is 



TABLE I. - SUMMARY OF CASES PRESENTED IN APPENDIX B FOR WHICH l 3 = 0 


Case Range of variables Nature of solution y(t) 



^0 

l 

1 

1 2 

b 


b 2 . 

4c 

0 

>> 

t 

>» 

111 

N 


A 

Arbitrary 3, 

Arbitrary 

> 0 

Arbitrary 

>0 

< 0 

> 0 

Periodic, positively infinite maxima 

B 


< 0 

Arbitrary 

> 0 

< 0 

< 0 

Periodic, negatively infinite minima 

C 


Arbitrary 


0 

3 

0 

> 0 

Monotone decreasing, unbounded 

D 


> 0 

> 0 

b 2 / 4 



> 0 

Periodic, positively infinite maxima 

E 


> 0 

< 0 




> 0 

Monotone increasing, bounded 

F 


< 0 

> 0 




< 0 

Monotone decreasing, bounded 

G 



< 0 

< 0 




< 0 

Periodic, negatively infinite minima 

H 




> 0 

> 0 

3 

> 0 

>-b 

Monotone, singular 

I 




< 0 

> 0 




Arbitrary 

Monotone, singular 

J 




> 0 

< 0 





Arbitrary 

Monotone, singular 

K 




< 0 

< 0 





<-b 

Monotone, singular 

L 




> 0 

Arbi 

trary 

Arbi 

trary 



z > o> z 3 > z 2 

Periodic, positively infinite maxima 

M 




< 0 







o> z > z 3 > z 2 

Periodic, bounded 

N 




> 0 







z 3 > Z 2 > z > 0 

Periodic, bounded 

0 




< 0 







Z 3 >Z 2 >0 - Z 

Periodic, negatively infinite minima 

P 




< 0 



, 




z 3 > z > o> z 2 

Periodic, bounded 

Q 




> 0 







Z 3 0 ^ Z Z 2 

Periodic, bounded 

R 




0 

— 




... 

Arbitrary 

Periodic solutions possible 

S 


0 

0 

— 



— 

... 

Arbitrary 

Periodic solutions not possible 


a The parameter so designated must, however, be restricted within a range appropriate to the limitations imposed on the other 
parameters for that case. 



(45) 


k 


2 


_ Z 2 _ b + (b 2 - 4c) 1 / 2 

Zg / 2 \l/2 

d b - (b - 4c) 


The solution given by equation (43) has a peak-to-peak period of u^(T) = 2K(k). The 
period of oscillation is then given by 


T = 


4K(k) 


Ai/2 




-b + (b 2 - 4c) 


1 / 21 1 / 2 \ l 2 J 


(46) 


Equation (43) does represent a solution of physical interest, since it is both periodic and 
finite. 


APPLICATION AND LIMITATIONS OF THE METHOD 

A requirement for application of the method presented in the preceding section is that 
Ax(t) « It is not easy to determine whether this inequality holds when a particular 
example of equations (1) and (2) is under consideration. If the equations correctly 
represent an actual physical situation and it is known that one of the two variables is 
oscillating by this mechanism with a small amplitude with respect to its mean value, this 
information may be used to justify application of the analytical procedure developed in 
this paper. This "physical" justification of the use of the small -amplitude approximation 
is appropriate, for example, to the equations describing the oscillations in a lightly 
ionized gas (ref. 11). 

Circumstances can arise, however, in which equations (1) and (2) will be given with 
a particular set of numerical coefficients and the nature of their solutions will not be 
known a priori. When this is the case, either the equations can be solved numer- 
ically on a computer to verify that Ax(t) « x-p or the "phase-plane” method of deter- 
mining the limits of x a &d y can be applied. The former procedure is exact but rend- 
ers a closed-form solution to the equations unnecessary for the particular numerical 
values of the coefficients Aj and C. so calculated. The phase -plane method presumes 
a knowledge of the literature on the class of equations represented by equations (1) and (2). 
This method is discussed in general terms in references 2 and 6. Application of this 
method to specialized cases of equations (1) and (2) (when certain of the A. and/or CL 
are zero) is discussed in the journal literature, of which references 1, 3, 5, 7, 8, and 
11 are examples. 
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The phase-plane method has been developed in the literature to determine the nature 
of the solutions to coupled sets of equations such as (1) and (2). This method pre- 
dicts whether the solutions are periodic, steady-state, growing, or damped. This 
method determines the trajectory that the variables x and y trace out on the 
X - y (phase) plane. The equation for this phase-plane trajectory is derived by 
dividing equation (2) by equation (1), which eliminates time as an explicit variable 
and yields 


dy _ A o + A ix + A 2 y + A 3 xy + A 4 x 2 +A 5 y 2 (4?) 

dx C 0 + C 1 x + C 2 y + C 3 xy + C 4 x 2 + C 5 y 2 

The requirement that the x(t) and y(t) be periodic requires the existence of at least one 
closed trajectory on the phase plane. If more than one exist, the possibility of stable 
oscillations of differing amplitudes is implied. The peak-to-peak amplitude of these 
curves in the x -direction is determined by the initial conditions and the zeros of the 
denominator; the peak-to-peak amplitude of y is determined by the initial conditions and 
the zeros of the numerator. 

The requirement that Ax(t) « x^ and that the solutions be periodic implies that 
the phase trajectories be narrow closed curves and that their width be small compared 
with the mean value of x> Xp at the center of the curve. This situation on the phase 
plane is illustrated schematically in figure 2. The closed-form solutions given in refer- 
ences 3 and 5 are appropriate only for small closed curves about x^ and yQ 4 , where the 
peak-to-peak amplitudes of x an d y are both small by comparison with their mean 
values. A further discussion of the phase-plane technique appears in references 2 
and 6. 

The peak-to-peak amplitudes of x(t) and y(t) are determined jointly by the initial 
conditions and the values of the coefficients A. and (L . For a given set of initial con- 
ditions on x and y, the requirement that Ax(t) « x^ places restrictions on the values 



Figure 2. - Schematic illustration of phase-plane 
trajectory. 
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of Aj and for which the present procedure will yield valid closed -form solutions. 
The restrictions on and implied by Ay(t) « Xj exclude many cases of interest. 
The advantage of the procedure developed herein is that the requirement of references 3 
and 5, that both Ax(t) « Xi and Ay(t) « y Q1 , is clearly more restrictive of the coeffi- 
cients A^ and than the present case, where only one amplitude is assumed to be small. 

If the coefficients of equations (1) and (2) are such that Ay(t) « Xj> the method of 
solution given herein may nonetheless fail if the time -dependent terms in the brackets of 
equation (6) are not small by comparison with the constant terms within these brackets. 
The latter terms are given by the l ^ in equations (10) to (13). Applying the procedure 
developed herein requires checking not only to see whether the coefficients are such that 
Ax(t) « Xp but also to see whether the following inequalities are satisfied: 

l 0 » Ax(t)(A 1 C 1 + AjA 2 + AqA 3 + 2A 4 C q )+ 2 x x Ax(t)(A x C 4 + A^ 

+ A X A 3 + 2A 4 Cj) + 3 X j Ax(t)(A 3 A 4 + 2A 4 C 4 ) (48) 

l 4 » Ax(t)^ jCg + 2A 2 A 3 + AgC j + 2A 4 C 2 + 2A & A^ 

+ 2 Xl A X (t)(A 3 C 4 + A§ + 2A 4 C 3 + 2A 4 A 5 ) (49) 

and 

l 2 » A X (t)(A 3 C 3 + 3A 3 A 5 + 2A 4 C 5 ) (50) 

where l Q , l j, and are given by equations (10), (11), and (12), respectively. The 
better these inequalities are satisfied, the more nearly will equation (6) approximate the 
Jacobian elliptic equation with constant coefficients. The solutions of equation (9) will 
then become better and better approximations to the small -amplitude solutions of equa- 
tions (1) and (2). 

This procedure has a further difficulty in common with all procedures involving a 
small -amplitude approximation. There is no assurance that oscillatory (or monotone) 
solutions to the small -amplitude approximation given herein will, for the same A C., 
and Xp correspond to oscillatory (or monotone) solutions to the exact equations (1) 
and (2). This state of affairs comes about because the time-dependent terms in the 
however small they may be, may possibly introduce damping into solutions whose small - 
amplitude approximation is periodic, or make periodic a solution whose small -amplitude 
approximation is monotone. 
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CONCLUSIONS 


A method has been developed which provides approximate solutions for the coupled 
set of first-order, nonlinear equations 

X = C Q + Cj x + c 2 y + c 3 xy + c 4 x 2 + c 5 y 2 (1) 

y = A Q + A x x + A 2 y + Ag xy + A 4 x 2 + Ag y 2 (2) 

where x and y are variables and Aj and are constants. The method provides 
solutions for y(t) and requires that the fluctuations in the other variable x(t) be small 
compared with its mean value. Additional requirements for the applicability of the 
method are set forth. This method is of particular value in obtaining oscillatory solutions 
for y(t), for which it yields amplitude as well as period and waveform. In general, the 

solutions for y(t) are given in terms of Jacobian elliptic functions. These functions may 

be approximately sinusoidal under certain conditions. The waveform and amplitude of 
X(t) are not given by this procedure. 

Classes of solutions are described for various combinations of the coefficients C. 
and Aj. The most general case, with Cg and Ag not equal to zero, is left in the form 
of an integral. With these two coefficients set equal to zero, nineteen different combina- 
tions of coefficients were treated. Among these combinations, five permitted finite - 
amplitude oscillatory solutions. 

The analysis presented herein extends the solutions available in the literature 
(refs. 3 and 5) to the case in which the peak-to-peak amplitude of y(t) may be compar- 
able to its mean value. In appendix C it is shown that the results obtained in this analysis 
reduce to those given in references 3 and 5 in the limit of small peak-to-peak fluctuations 
in y(t). The present analysis therefore contains the previously obtained closed-form 
solution to equations (1) and (2) as a special case. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, November 14, 1967, 

129-02-03-05-22. 
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APPENDIX A 


A i 


b 

C, 


CN 


c 

f 

g 

h 

i 

j 

K(k) 

k 


k’ 


^ 3 
SN 


SYMBOLS 


constant parameters appearing 
in eq. (2) 

constant defined by eq. (33) 

constant parameters appearing 
in eq. (1) 

Jacobian elliptic cosine function 
(see ref. 13) 

constant defined by eq. (34) 
constant defined by eq. (24) 
constant parameter 
constant defined by eq. (25) 

constant defined by eq. (23) 

complete elliptic integral of 
first kind 

modulus of Jacobian elliptic 
functions 

complementary elliptic modulus, 
o 

equal to 1 - k 

constant defined by eq. (10) 

constant defined by eq. (11) 

constant defined by eq. (12) 

constant defined by eq. (13) 

Jacobian elliptic sine function 
(see ref. 13) 


T 

TN 

t 

U 1 


y 

^max 

y o 

y 00 

y 01 

z 

Z max 

6 


X 

AX 

*1 

*2 

O) 


peak -to -peak period of oscilla- 
tion, sec 

Jacobian elliptic tangent (see 
ref. 13) 

time, sec 

argument of elliptic function, 
incomplete elliptic integral of 
first kind 

variable appearing in eq. (2) 

maximum value of y 

minimum value of y 

mean value of y 

amplitude of fluctuation in y 

relative amplitude of y, defined 
by eq. (16) 

peak -to -peak amplitude of y, 
defined by eq. (19) 

small perturbation amplitude, 
defined by eq. (3) 

frequency, Hz 

variable appearing in eq. (1) 
peak -to -peak amplitude of x 
mean value of x 
amplitude of fluctuation in x 
frequency, rad/sec 
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APPENDIX B 


SOLUTIONS FOR y(t) WHEN l 3 = 0 

The algebraic manipulations required to obtain closed-form solutions to equation (9) 
can be somewhat tedious and may involve mathematical functions unfamiliar to the reader 
It was therefore deemed desirable to exhibit the solutions for the case £3 = 0 and, thus 
make it easier for the reader to obtain y(t) from a knowledge of the coefficients Aj and 
Cj, or vice-versa. The solutions for the cases listed in table I (p. 17) are therefore 
presented here. 


Case l 2 t 0 


Case A, l ^ > 0, (b^ - 4c) <0, Z > 0: 

9 

When l g is greater than zero and the discriminant b - 4c is less than zero, 
the roots Zg and Zg are complex conjugates, and the lower limit in the range 
of integration is equal to the only real root of the denominator, Zj = 0. Under 
these circumstances, the integral in equation (40) is given by formula 239:00 of 
reference 13, which is 


Ui(k, t) = - 

g 



where 


g 


= c 


-1/4 


the amplitude Z has the functional form 


- 1 / 2 , 


Z = c 


fl - CNu, 


^1 + CNUj 


| = y -y 0 


2 

And the elliptic modulus k has the value 



(41) 

(Bl) 


(B2) 


(B3) 


(B4) 
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o 

The requirement that k be real and positive places the following requirements on b 
and c 


c > 0 


(B5) 


• 1 < 


2c 


1/2 


< 1 


(B6) 


An examination of equation (B3) reveals that Z has the period u^(T) = 4K(k), where K 
is the complete elliptic integral of the first kind. Therefore, 

T , 4K(k) / 3 
c V4 \»2 

However, during one period the amplitude will become infinite when the elliptic argument 
u i is some multiple of 2K(k). Thus, while the solutions given by equation (B3) are 
periodic, they involve infinite amplitudes, which are unrealistic for the solution to a 
physical problem and beyond the scope of the small -amplitude approximation. 

Case B, l 2 <0, (b 2 - 4c) <0, Z < 0: 

This case is identical to the previous one, except that the argument of CN is imagin- 
ary and Z is negative. 


, 1/2 


sec 


(B7) 


u 2 (k',t) = iujfc^t) 



(B8) 


In reference 13, expression 125. 02 is a formula for the CN function of imaginary argu- 
ment, where 


CN i(uj(k)) = [CNU^')]- 1 


and 


k' 2 = 1 - k 2 



(B9) 


(BIO) 
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If Uj is now understood to be a function of k’, equation (B7) for the amplitude of Z may 
be written as 


Z = c 


1/2 


yCNUj + 1 J 


-yo 


(Bll) 


Thus, if 1 2 < 0, the amplitude Z is infinite when = 2K(k'). Equation (Bll) may also 
be obtained from formula 243. 00 in reference 13. Therefore, case B also does not lead 
to physically realistic solutions within the scope of the small -amplitude approximation. 

o 

Case C, l g arbitrary, b - 4c = 0, b = c = 0: 

In this case, the integral of equation (40) assumes the much simpler form 


21 , 


1/2 


/ 

**0 


dZ 

,3/2 


(B12) 


This integral is obviously monotone and divergent, regardless of the value of l ^ CaseC 
therefore does not result in bounded periodic solutions. 

Case D, 1 2 > 0, b 2 - 4c = 0, b < 0, Z > 0: 

In this case Zg = Zg = -b/2 and the basic integral (eq. (40)) is given by 


( 21 , 


1/2 


dZ 


,3/2 


I z i/2 (z*D vw 


-i 

b/ 


tan -i 


(B13) 


The amplitude of the function Z is therefore 


Z = b tan 2 


ll 2 b\ 

\u) 


1/2- 


y -y 0 


(B14) 


The function Z will have a peak-to-peak period of i r radians, which results in a period of 
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The amplitude of Z becomes infinite, however, and so equation (B14) cannot represent 
physically realistic solutions within the scope of the small amplitude approximation. 

Case E, 1 2 > 0, (b 2 - 4c) = 0, b < 0, Z > 0: 

When the constant b is negative, the integral in equation (B13) is given by 



This solution is bounded, monotone, and a physically realistic solution to the problem. 

Case F, 1 2 < °, (b 2 - 4c) = 0, b > 0, Z < 0: 

When the identity tan iy = i tanh x is used, the solution for this case is obtained 
from equation (B14) as follows: 



(B17) 


This solution is not of interest in the present analysis since it is nonperiodic. 

Case G, l 2 < 0, (b 2 - 4c) - 0, b < 0, Z < 0: 

The functional form for this case can be obtained by using the identity tanh ix = i tan x 
in equation (B16) to yield 



(B18) 


This solution has the period given by equation (B15), but it is not physically realistic 
since it has infinite amplitudes and is therefore beyond the scope of the small -amplitude 
approximation. 
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Case H, Z 2 > 0, (b 2 - 4c) > 0, c = 0, b > 0: 

From equations (35) and (40) it is apparent that in this case the basic integral is 
given by 




(Z +b) 1 / 2 -b 1 / 2 Z 
(Z + b) 1//2 + b 1 / 2 0 


(B19) 


Because a singularity exists at the lower limit of integration, this case cannot represent 
a small -amplitude solution to the differential equation. 

Case I, Z 2 < °, (b 2 - 4c) > 0, c = 0, b > 0: 

The negative value of l 2 results in the left side of equation (B19) being multiplied by 
V-i- This process does not remove the singularity at Z = 0; therefore, this case also 
cannot represent a small -amplitude solution. 

Case J, Z 2 > 0, (b 2 - 4c) > 0, c = 0, b < 0: 
and 

Case K, Z 2 < °, (b 2 - 4c) > 0, c = 0, b < 0: 

In both these cases, the basic integral, given by equation (B19), is divergent at its 
lower limit. Therefore, the basic differential equation does not have a small -amplitude 
solution. 

Case L, Z g > 0, (b 2 - 4c) > 0, Z > 0 > Zg > Z 2 : 

In this case, it is assumed that the parameter Z 9 is positive, that the discriminant 

2 “ 
b - 4c is positive, and that the constant parameters are such that the roots Zg and 

Z 2 , given by equations (38) and (39), are both negative. In this case, the value of the 

basic integral (eq. (40)) is given by formula 237. 00 of reference 13. 

The solution in this case is 


uj(k,t) = - 
g 



(Bl) 
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where 


,3/2 



The amplitude of Z has the functional form 

/ 2 X 1 / 2 

Z = b ~ \ b ~ 4c ) . TN 2 Uj(k) = y - y Q 

2 

and the modulus of the elliptic function is equal to 


(B20) 


(B21) 


k 


2 



(B22) 


The requirement that 0 < k 2 < 1 is automatically satisfied when both roots Zg and Zg 
are negative. Equation (B21) has a peak-to-peak period of u^T) = 2K(k), which results 
in the following period of oscillation: 


T = 


b + b 


4K(k) 

\l/2 

_ An I 

) 


ll/2 



It is clear from equation (B21) that Z becomes infinite at = K, 3K, etc. 
this solution cannot represent a small -amplitude solution. 


(B23) 


Therefore, 


Case M, 1 2 < 0, (b 2 - 4c) > 0, 0 > Z > Zg > Z 2 : 

This case is closely related to the previous one, in which equation (Bl) implies that 


u 2 (k,t) = iu 1 (k,t) 


(B24) 


Using formula 125. 02 of reference 13 yields 


28 


(B25) 


TNiu 1 (k) = iSNUjfc') 


where the complementary modulus is given by 


k ,2 _ b - (b 2 - 4c) 


1/2 


b + (b 2 - 4c) 

and equation (B21) becomes for this case 


\l/2 


Z = - 


(b 2 - 4c) 


1/2 


SN^u 1 (k’) =y - y Q 


(B26) 


(B27) 


The period is given by equation (B23) with k replaced by k' . The amplitude Z is finite, 
negative, and periodic; therefore, equation (B27) is a physically interesting solution to 
the problem. Equation (B27) could also have been obtained from expression 236. 00 of 
reference 13. 

Case N, >0, (b 2 - 4c) >0, Zg > Z g > Z > 0: 

When both of the roots given by equations (38) and (39) are positive and the range of 
Z is confined to Zg ^ Z < 0, the solution to the basic integral of equation (40) is given 
by formula 233. 00 of reference 13. 


Ui(k,t) = - 

g 



where 


g = 


,3/2 


- b + ^b 2 - 4c) 


1/2 


n 1/2 


The amplitude of Z has the functional dependence 


Z = Z 2 SN 2 Ul (k) = - - b + (b 2 - 4c) 


\l/2 


SN u x (k) = y -y Q 


(Bl) 


(B28) 


(B29) 
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The total peak-to-peak amplitude of Z is then 


^max ^2 



(b 2 - 4c) 



and the modulus of the elliptic function is 


k 


2 



/ 2 N 1 / 2 

_ b + (b* - 4c ) 

(2 X 1 / 2 

b - (b^ - 4c) 


(B30) 


(B31) 


The solution given by equation (B29) has a peak-to-peak period of Uj(T) = 2K(k). The 
period of oscillation is then given by 


T = 


4K(k) 




/ 2 \ V 2 

- b + (b* - 4c) 


1 / 21 1 / 2 \ l 2 / 


(B32) 


Equation (B29) does represent a solution of physical interest, since it is both periodic 
and finite. 

Case O, l 2 < 0, (b 2 - 4c) >0, Zg > > 0 > Z; 

This case is identical to the previous one, except that the argument Uj is imaginary 
and Z < 0. Using formula 125. 02 in reference 13 yields 


SNiUj(k) = i TNu^(k') 


(B33) 


Equation (B29) may then be written 


Z = - Z 2 TN 2 Ul (k’) =- 


/ 

2 x 1 / 2 ! 

b + (fc 

r - 4c) 


where the complementary modulus is 


k’ 2 = 1 - k 2 


_ -2(b 2 - 4c) 1/2 
b - (b 2 - 4c) 


(B34) 


(B35) 
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I 


and the period is given by equation (B32) with k replaced by k' . Equation (B34) cannot 
represent a small -amplitude solution to the problem since the amplitude becomes infinite. 
Equation (B34) may also be obtained from expression 232. 00 of reference 13. 

Case P, Z g < 0, (b 2 - 4c) > 0, Zg > Z > 0 > Z 2 : 

In the case in which there is only one positive root and the range of integration is 
restricted between Zg > Z > 0, the solution to the basic integral of equation (40) is given 
by formula 235. 00 of reference 13. For this case, 


2 1 , 


Ul ( M ) =- .. _ 

g V 3 


(B36) 


where 


g = 


^b 2 - 4cj 

The amplitude of Z has the functional dependence 


1/4 


(B37) 


Z = 


Z 2 ZgSN 2 u 1 (k) 


-cSN u x (k) 


(„* - 4c) 1/2 - 1 

/ 

2 \ 1/2 

- b + b 

1 - 4c) 

\ / 2 

\ 

/ 


SN u x (k) 


(B38) 


where the modulus k is equal to 


k 2 = 


Zg-Z 2 2 


1 -■ 


^b 2 - 4c 


y/2 


(B39) 


The requirement that 0 < k“ < 1. 0 will be met, since c < 0 in this case. The total 
peak -to -peak amplitude of the fluctuations in Z is given by 


^max ^3 2 


- b + (b 2 - 4cj 


1/2 


(B40) 
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Equation (B38) has a peak-to-peak period of u^(T) + 2K(k), so that the period of the 
oscillation is given by 


T = 2K(k)g 




The solution of equation (B38) is positive, periodic, and finite. 


Case Q, 1 2 > 0, (b 2 - 4c) > 0, Zg > 0 > Z > Z 2 : 


(B41) 


The solution to this case is identical to the previous one, except that the argument 
u^(k) is imaginary. With the use of equation (B33), equation (B38) may be written as 


Z = 


-Z^gT^UjOt*) 

Z 2- Z 3- Z 3 TN S< k '> 


(B42) 


This solution is periodic and finite and is therefore a desired solution to the problem. 

The period of Z is given by equation (B41), with k replaced by k\ The complementary 
modulus k’ is equal to 1 - k , where k is given by equation (B39) and Z max is equal 
to Zg. Equation (B42) can also be obtained from expression 234. 00 of reference 13. 

Case R, l 2 = l^ = 0, 1 0: 

The next situation, of lesser mathematical complexity, arises when l ^ * 0, 
l 2 = Zg = 0, and l q may have any value. For this case, equation (9) may be written 

y = *o + *i y ( B43 ) 


In terms of the variable Z 


defined in equation (16), this equation may be written 


Z =Z 0 + *l y 0 + *l Z 


(B44) 


If 


2 




is defined as 


2 


Cl) 



(B45) 
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equation (B44) may be written as 


Z + 


u> 


' z =*o- 


w 2 y 0 


(B46) 


The solution to this equation is 


Z = 



(cos cot - 1) 


(B47) 


which satisfies the boundary conditions of equation (18). The requirement that the solu- 
tions be periodic implies that 


and therefore, 


< 0 


co 2 > 0 


(B48) 


And the requirement that the amplitude of y(t) be either zero or positive places the 
following restriction on the permissible value of Iq with respect to ^ and Yq: 

= 2 y 0 w2 (B49) 


The constant l ^ must therefore be positive and larger than indicated in equation (B49) if 
periodic motion which satisfies the boundary conditions and has a positive amplitude is to 
occur. 

The period of the fluctuation in y is 


T = — = — (B50) 


and the total amplitude of the fluctuation in y is 

l 


Av = Z =2 
y max 


v - 0 

70 ~~2 
Ct> 


= 2 


0 

y ° + r 


(B51) 
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Case S, ^ = ^ 2=^3 = °: 

The most elementary special case occurs when only Iq is nonzero. In this case, 
equation (9) may be written 


y = z = l o 

Integrating twice and applying the boimdary conditions of equation (18) gives 

z = iot 2 

2 


(B52) 


(B53) 


There are no periodic solutions to equation (B52), and the only solution is monotone in 
time. 
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APPENDIX C 


COMPARISON WITH PRIOR LITERATURE 

One of the solutions to equations (1) and (2) which appears in appendix B can be com- 
pared with the solution of a much more restricted case that has appeared in the literature. 
Lotka (ref. 3) and Volterra (ref. 5) both considered periodic, small -perturbation solu- 
tions to the system of equations 


x = c x x + c 3 xy 

(Cl) 

y = A 2 y + A 3 xy 

(C2) 


where these equations have been written in the notation used herein. These equations 
differ from equations (1) and (2) in that 


and 



<C3) 


(C4) 


In the present analysis, it is assumed that the fluctuations in x are small by com- 
parison with the mean value of x, but that the peak-to-peak amplitude of y can be either 
smaller or larger than the mean value of y. Lotka and Volterra assumed that the per- 
turbations of both x and y were small with respect to their mean values. These 
authors expanded x and y in the form 


X = Xi + X 2 (t) 

(C5a) 

Xi » X 2 (t) 

(C5b) 

y = yoo + yoi (t) 

(C6a) 

yoo »yoi (t) 

(C6b) 
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where Xj and y Q0 , the time -independent mean values of x and y, are defined as 


= - -f- (C7) 

A 3 

y 00 * - ^ <C8) 

Substituting equations (C5a), (C6a), (C7), and (C8) into equations (Cl) and (C2) yields 

x 2 =C 3 y 0l (x l + x 2 } (C9) 

y 01 =A 3 x 2 (y 00 + y 01 ) ( C1 °) 

Differentiating equation (CIO) with respect to time and eliminating the first order time 
derivatives yield the second-order equation 


y 01 _A 3 C 3 y 01^ x l + x 2^ y 00 + y 01^ + A 3 X 2^00 + y 01^ ( C11 ) 

At this point the approximations given by equations (C5b) and (C6b) must be intro- 
duced. When these approximations are introduced and the products of small quantities 
are neglected, equation (Cll) becomes 


y 01 A 3 C 3 x l y 00 y 01 


(C12) 


By making use of equations (C7) and (C8), the frequency may be written 

w = (' A 3 C 3 x l y Oo) 1/2 = (~ A 2 C l) 1/2 ( C13 ) 

A solution of equation (C12) which is consistent with the initial conditions of equation (4) 
is 


y Q i(t) = 6 cos cut 


(C14) 


where 6 is (in this appendix) the time -independent amplitude of the perturbation in y. 
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The waveform of the oscillation in y is then obtained by substituting equation (C14) into 
equation (C6a), which gives 


y(t) = y 0 o + 6 cos wt 


(C 15) 


where 


y o = y oo + 6 


(C16) 


is considered to be the maximum value of y and to occur at t = 0. The frequency of the 
oscillations is, from equation (C13), 

(- A 3 C 3^1 y Oo) 1/2 (" A 2 C l) 1/2 

v = V 3 3 1 jgg = V 2 ~.V. cps (Hz) (C 17) 

2ir 2ir 

and the period is 

T = sec (C18) 

( _A 3 C 3 x l y Oo) 1//2 ( _A 2 C l) 1//2 

which is the result obtained by Lotka (ref. 3) and Volterra (ref. 5). 

This prior result will be compared with the results of the present analysis, in the 
limiting situation given by equations (C5b) and (C6b). Substituting equations (C3)and (C4) 
into equations (10) to (13) yields 


J 0 = Z 3 = 0 (C19) 

1 1 = ( A 2 + A 3 x l) 2 + A 3 C 1 X 1 (C20) 

Z 2 =A 3 C 3 x 1 ( C21 ) 


If x(0) = Xj is chosen as an initial condition, the initial conditions on y(t) given by 
equation (4) then imply 


A 2+ A 3 x i = ° (C22) 

so that equation (C20) becomes 
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*1 =A 3 C 1*1 


(C23) 


Since 1 3 = 0 , the results presented in appendix B are appropriate. In actual physical 
applications, x and y will represent the magnitude of real physical quantities and will 
be positive. As can be seen by manipulating equations (Cl) and (C2), either Ag or Cg 
must be negative in order for periodic solutions to exist. The parameter in 
equation (C21) will therefore be negative. As in the analysis leading to equation (C15), it 
is assumed here that yQ is the maximum value of y, occurring at t = 0. This assump- 
tion will result in the ancillary variable Z, defined by equation (16), being negative. The 
three conditions, l g = 0, l^< 0, and Z < 0, therefore imply that case M of appendix B 
is appropriate for comparison with the results of Lotka and Volterra. If a positive par- 
ameter 


V ~~ 


c 3yo 


(C24) 


is defined, the parameters b and c given by equations (33) and (34) are equal to 


b 2 y 0 (2 “ ^ 


(C25) 


c = 3y q (1 - 77 ) 


(C26) 


Equations (38) and (39) for the roots are 


r 


z 2 = i^ 


(2 - v) + 


H Y v+2) _ 


1/2 


r 


^3 = ~ 4 0’ 


(2 - u) - 


(” ■ I)' 1 ' + 2> 


1/2 


(C27) 


(C 28 ) 


From equation (B27), it can be seen that the waveform predicted by the present analysis 
is 


y(t) ~ y Q + ZgSlN a 1 (k f ) 


(C29) 
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An inspection of equation (C28) reveals that a necessary and sufficient condition for the 
small -amplitude limit Zg « y Q is 77 - 1. 0. As 77 approaches 1 as a limit, equa- 
tions (C24) to (C28) become 


C 3 y 0 » - c x 

(C30) 

b “5 y ° 

(C31) 

c « 0 

(C32) 

Z 2“-^ y 0 

(C33) 

z 3 « 0 

(C34) 


In this limit, the period of oscillation given by equation (B23) is equal to 


T ~ 4^0*’) = 4K(k’) 

~ M ' 2 (-w s ^) 1/2 


From equation (B26), in the limit 77 — 1.0, 



(C35) 


(C36) 


Thus, when 77 — 1. 0, 

K(k’) « 1 (C37) 

2 

Substituting equations (C37), (C30), and (C22) into equation (C35) shows the period of 
oscillation to be 


T « 


2jt 

(- y 0 Xi A 3 C 3 ) l/2 



(C38) 
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which, if Yq « Yqq, is in agreement with the period derived by Lotka and Volterra in 
equation (C18). 

2 

As the elliptic modulus k' approaches zero, the Jacobian elliptic sine becomes 
equal to the conventional trigonometric sine (ref. 13). 

SNu 1 (k*) - sinu 1 (0) (C39) 

The argument u^(k' — 0) may be obtained from equations (B19) and (B20). 

u i<°) “ j t (-y 0 x 1 A 3 c 3 ) 1/2 = s£ (c40) 

When this argument is substituted into equation (C29) the small -amplitude limit is 

z 

y(t) = y Q + Z 3 sin 2 — = y 0 + — (1 - cos out) (C41) 

2 2 

Identifying Zg with the small amplitude 5, 

Z 3 = -26 (C42) 

and using equation (C16) to rewrite equation (C41) yield the limiting small -amplitude 
expression 


y(t) « y Q0 + 6 cos o>t (C43) 

Thus, the frequency and waveform predicted by this analysis in the small -amplitude 
variation of both x(t) and y(t) agree with the small -perturbation results of Lotka and 
Volterra, which are given by equation (C15). 
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